//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testFissionSource.cc
 * \author Jeremy Roberts
 * \date   Sep 4, 2011
 * \brief  Test of class FissionSource.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//


#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "test/Testing.hh"
#include "test/TestInput.hh"

#include "angle/Definitions.hh"
#include "angle/MomentsIndex.hh"
#include "angle/AngleTraits.hh"

#include "material/BuiltInCrossSections.hh"
#include "material/Materials.hh"

#include "tran/Moments_Field.hh"
#include "tran/CartMesh.hh"
#include "tran/Traits.hh"
#include "tran/Definitions.hh"
#include "tran/FissionSource.hh"

#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/SP.hh"

// Equations
#include "tran/DD_1D_Equations.hh"
#include "tran/DD_2D_Equations.hh"
#include "tran/DD_3D_Equations.hh"

using namespace std;
using namespace slabtran;
using namespace util;

// typedefs
typedef FissionSource<DD_1D_Equations,_1D>::SP_fissionsource SP_source_1d;
typedef FissionSource<DD_2D_Equations,_2D>::SP_fissionsource SP_source_2d;
typedef FissionSource<DD_3D_Equations,_3D>::SP_fissionsource SP_source_3d;

// 1d tests.
int test1d_construct(SP_source_1d &source);
int test1d_build_fission_density(SP_source_1d &source);
int test1d_build_fission_source(SP_source_1d &source);
int test1d_verify_source(SP_source_1d &source);
// 2d tests
int test2d_construct(SP_source_2d &source);
int test2d_build_fission_density(SP_source_2d &source);
int test2d_build_fission_source(SP_source_2d &source);
int test2d_verify_source(SP_source_2d &source);
// 3d tests
int test3d_construct(SP_source_3d &source);
int test3d_build_fission_density(SP_source_3d &source);
int test3d_build_fission_source(SP_source_3d &source);
int test3d_verify_source(SP_source_3d &source);

int main(int argc, char *argv[])
{
  Require(argc==2);
  int test = atoi(argv[1]);
  cout << " ------ testFissionSource ------ " << endl;
  int test_value = 0;
  int number_of_tests = 13;

  // Declare sweep sources for tests.
  SP_source_1d s1; // 1-d fission source
  SP_source_2d s2; // 2-d fission source
  SP_source_3d s3; // 3-d fission source

  if ((test >=  1 && test <=  4) || test == number_of_tests) test_value += test1d_construct(s1);
  if ((test >=  2 && test <=  4) || test == number_of_tests) test_value += test1d_build_fission_density(s1);
  if (test ==  3 || test ==  4   || test == number_of_tests) test_value += test1d_build_fission_source(s1);
  if (test ==  4 || test == number_of_tests                ) test_value += test1d_verify_source(s1);

  if ((test >=  5 && test <=  8) || test == number_of_tests) test_value += test2d_construct(s2);
  if ((test >=  6 && test <=  8) || test == number_of_tests) test_value += test2d_build_fission_density(s2);
  if (test ==  7 || test ==  8   || test == number_of_tests) test_value += test2d_build_fission_source(s2);
  if (test ==  8 || test == number_of_tests                ) test_value += test2d_verify_source(s2);

  if ((test >=  9 && test <= 12) || test == number_of_tests) test_value += test3d_construct(s3);
  if ((test >= 10 && test <= 12) || test == number_of_tests) test_value += test3d_build_fission_density(s3);
  if (test == 11 || test == 12   || test == number_of_tests) test_value += test3d_build_fission_source(s3);
  if (test == 12 || test == number_of_tests                ) test_value += test3d_verify_source(s3);

  if (test > number_of_tests || test < 1)
    cout << " test " << test << " not available. " << endl;
  if (test_value)
    cout << " failed with value = " << test_value << endl;
  else
    cout << " passed." << endl;
  return test_value;
}

//===========================================================================//
// Reference data
//===========================================================================//
// Set the reference fission density.
double f_ref[] = {
// material 0
           0.00000 * 0.00000 * 1.00  +  // SigmaF_1 * nu_1 * phi_1
           0.00000 * 0.00000 * 1.00  ,  // SigmaF_2 * nu_2 * phi_2
// material 1
           0.00268 * 2.50000 * 1.00  +  // SigmaF_1 * nu_1 * phi_1
           0.04964 * 2.50000 * 1.00  ,  // SigmaF_2 * nu_2 * phi_2
// material 2
           0.00312 * 2.50000 * 1.00  +  // SigmaF_1 * nu_1 * phi_1
           0.06168 * 2.50000 * 1.00  ,  // SigmaF_2 * nu_2 * phi_2
// material 3
           0.00224 * 2.50000 * 1.00  +  // SigmaF_1 * nu_1 * phi_1
           0.00748 * 2.50000 * 1.00  }; // SigmaF_2 * nu_2 * phi_2

//===========================================================================//
// 1-D TESTS
//===========================================================================//
int test1d_construct(SP_source_1d &source)
{
  // Make test input.
  TestInput::SP_input input;
  input = new TestInput(2,          // number groups
                        0,          // legendre order
                        1,          // dimension
                        4,          // quad order
                        GAUSSLEG);  // quad type

  // Make a simple 1d mesh
  Mesh<_1D>::SP_mesh mesh;
  mesh = new Mesh<_1D>(10.0,        // width [cm]
                       8);          // meshes
  // Set material id's.
  Mesh<_1D>::Vec_Id matids(mesh->number_cells(), 1);   // fuel I
  matids[0] = 0;                                       // water
  matids[1] = 2;                                       // fuel II
  matids[2] = 3;                                       // fuel II+Gd
  mesh->set_cell_id(MATERIAL_ID, matids);

  // Make a simple angular mesh.
  AngularMesh<_1D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_1D>(input->quadrature_order(),
                                     input->quadrature_type());

  // Make a state and set the moments to unity.
  State<DD_1D_Equations,_1D>::SP_state state;
  state = new State<DD_1D_Equations,_1D>(input, mesh, angularmesh);
  for (int g = 0; g < state->number_groups(); g++)
  {
    State<DD_1D_Equations,_1D>::Moments_Field_t &phi = state->get_phi(g);
    for (int i = 0; i < phi.size(); i++)
      phi[i] = 1.0;
  }

  // Use a built-in library.
  BuiltInCrossSections xs(BuiltInCrossSections::XS2G4M);
  Materials::SP_materials materials;
  materials = xs.build();

  // Build the fission source.
  source = new FissionSource<DD_1D_Equations,_1D>(state,
                                                  mesh,
                                                  materials);

  cout << " ... passed 1-d FissionSource construct test. " << endl;
  return 0;
}

int test1d_build_fission_density(SP_source_1d &source)
{
  Require(source);

  // Build the fission density for keff = 1.1
  source->build_fission_density(1.1);

  // Get the density.
  const FissionSource<DD_1D_Equations,_1D>::FissionDensityField &f =
      source->get_fission_density();
  FissionSource<DD_1D_Equations,_1D>::SP_mesh mesh = source->mesh();
  // Verify length of density vector
  TEST(f.size() == mesh->number_cells());

  // Check fission density.
  const Mesh<_1D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(f(i, 0), f_ref[matid[i]]));

  cout << " ... passed 1-d FissionSource build_fission_density test. " << endl;
  return 0;
}

int test1d_build_fission_source(SP_source_1d &source)
{
  Require(source);

  // Build the fission density for group 1 (indexed as 0!!)
  source->build_fission_source(0);

  // Get the source.
  const FissionSource<DD_1D_Equations,_1D>::FissionSourceField &f =
      source->get_fission_source();
  FissionSource<DD_1D_Equations,_1D>::SP_mesh mesh = source->mesh();
  // Verify length of density vector
  TEST(f.size() == mesh->number_cells());

  // Set reference factor.  Chi is 1, so factor is just 1/(2*keff) for 1d.
  double factor = 0.5 / 1.1;

  // Check factor.
  TEST(source->factor() == factor);

  // Check fission source vector directly.
  const Mesh<_1D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(f(i, 0), factor * f_ref[matid[i]]));

  cout << " ... passed 1-d FissionSource build_fission_source test. " << endl;
  return 0;
}

int test1d_verify_source(SP_source_1d &source)
{
  Require(source);

  // Build the fission density for group 1 (indexed as 0!!)
  source->build_fission_source(0);

  // Get the mesh.
  FissionSource<DD_1D_Equations,_1D>::SP_mesh mesh = source->mesh();

  // Set reference factor.  Chi is 1, so factor is just 1/(2*keff) for 1d.
  double factor = 0.5 / 1.1;

  // Check factor.
  TEST(source->factor() == factor);

  // Check fission source using q_f interface.
  const Mesh<_1D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(source->q_f(i, 0, 0), factor * f_ref[matid[i]]));

  cout << " ... passed 1-d FissionSource verify_source test. " << endl;
  return 0;
}

//===========================================================================//
// 2-D TESTS
//===========================================================================//
int test2d_construct(SP_source_2d &source)
{
  // Make test input.
  TestInput::SP_input input;
  input = new TestInput(2,          // number groups
                        0,          // legendre order
                        2,          // dimension
                        4,          // quad order
                        LEVELSYM);  // quad type

  // Make a simple 2d mesh
  Mesh<_2D>::SP_mesh mesh;
  mesh = new Mesh<_2D>(10.0,        // width [cm]
                       8);          // meshes
  // Set material id's.
  Mesh<_2D>::Vec_Id matids(mesh->number_cells(), 1);   // fuel I
  matids[0] = 0;                                       // water
  matids[1] = 2;                                       // fuel II
  matids[2] = 3;                                       // fuel II+Gd
  mesh->set_cell_id(MATERIAL_ID, matids);

  // Make a simple angular mesh.
  AngularMesh<_2D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_2D>(input->quadrature_order(),
                                     input->quadrature_type());

  // Make a state and set the moments to unity.
  State<DD_2D_Equations,_2D>::SP_state state;
  state = new State<DD_2D_Equations,_2D>(input, mesh, angularmesh);
  for (int g = 0; g < state->number_groups(); g++)
  {
    State<DD_2D_Equations,_2D>::Moments_Field_t &phi = state->get_phi(g);
    for (int i = 0; i < phi.size(); i++)
      phi[i] = 1.0;
  }

  // Use a built-in library.
  BuiltInCrossSections xs(BuiltInCrossSections::XS2G4M);
  Materials::SP_materials materials;
  materials = xs.build();

  // Build the fission source.
  source = new FissionSource<DD_2D_Equations,_2D>(state,
                                                  mesh,
                                                  materials);

  cout << " ... passed 2-d FissionSource construct test. " << endl;
  return 0;
}

int test2d_build_fission_density(SP_source_2d &source)
{
  Require(source);

  // Build the fission density for keff = 1.1
  source->build_fission_density(1.1);

  // Get the density.
  const FissionSource<DD_2D_Equations,_2D>::FissionDensityField &f =
      source->get_fission_density();
  FissionSource<DD_2D_Equations,_2D>::SP_mesh mesh = source->mesh();
  // Verify length of density vector
  TEST(f.size() == mesh->number_cells());

  // Check fission density.
  const Mesh<_2D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(f[i], f_ref[matid[i]]));

  cout << " ... passed 2-d FissionSource build_fission_density test. " << endl;
  return 0;
}

int test2d_build_fission_source(SP_source_2d &source)
{
  Require(source);

  // Build the fission density for group 1 (indexed as 0!!)
  source->build_fission_source(0);

  // Get the source.
  const FissionSource<DD_2D_Equations,_2D>::FissionSourceField &f =
      source->get_fission_source();
  FissionSource<DD_2D_Equations,_2D>::SP_mesh mesh = source->mesh();
  // Verify length of density vector
  TEST(f.size() == mesh->number_cells());

  // Set reference factor.  Chi is 1, so factor is just 1/(4pi*keff) for 2d.
  double factor = inv_four_pi / 1.1;

  // Check factor.
  TEST(source->factor() == factor);

  // Check fission source vector directly.
  const Mesh<_2D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(f[i], factor * f_ref[matid[i]]));

  cout << " ... passed 2-d FissionSource build_fission_source test. " << endl;
  return 0;
}

int test2d_verify_source(SP_source_2d &source)
{
  Require(source);

  // Build the fission density for group 1 (indexed as 0!!)
  source->build_fission_source(0);

  // Get the mesh.
  FissionSource<DD_2D_Equations,_2D>::SP_mesh mesh = source->mesh();

  // Set reference factor.  Chi is 1, so factor is just 1/(4pi*keff) for 2d.
  double factor = inv_four_pi / 1.1;

  // Check factor.
  TEST(source->factor() == factor);

  // Check fission source using q_f interface.
  const Mesh<_2D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(source->q_f(i, 0, 0), factor * f_ref[matid[i]]));

  cout << " ... passed 2-d FissionSource verify_source test. " << endl;
  return 0;
}

//===========================================================================//
// 3-D TESTS
//===========================================================================//
int test3d_construct(SP_source_3d &source)
{
  // Make test input.
  TestInput::SP_input input;
  input = new TestInput(2,          // number groups
                        0,          // legendre order
                        3,          // dimension
                        4,          // quad order
                        LEVELSYM);  // quad type

  // Make a simple 3d mesh
  Mesh<_3D>::SP_mesh mesh;
  mesh = new Mesh<_3D>(10.0,        // width [cm]
                       8);          // meshes
  // Set material id's.
  Mesh<_3D>::Vec_Id matids(mesh->number_cells(), 1);   // fuel I
  matids[0] = 0;                                       // water
  matids[1] = 2;                                       // fuel II
  matids[2] = 3;                                       // fuel II+Gd
  mesh->set_cell_id(MATERIAL_ID, matids);

  // Make a simple angular mesh.
  AngularMesh<_3D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_3D>(input->quadrature_order(),
                                     input->quadrature_type());

  // Make a state and set the moments to unity.
  State<DD_3D_Equations,_3D>::SP_state state;
  state = new State<DD_3D_Equations,_3D>(input, mesh, angularmesh);
  for (int g = 0; g < state->number_groups(); g++)
  {
    State<DD_3D_Equations,_3D>::Moments_Field_t &phi = state->get_phi(g);
    for (int i = 0; i < phi.size(); i++)
      phi[i] = 1.0;
  }

  // Use a built-in library.
  BuiltInCrossSections xs(BuiltInCrossSections::XS2G4M);
  Materials::SP_materials materials;
  materials = xs.build();

  // Build the fission source.
  source = new FissionSource<DD_3D_Equations,_3D>(state,
                                                  mesh,
                                                  materials);

  cout << " ... passed 3-d FissionSource construct test. " << endl;
  return 0;
}

int test3d_build_fission_density(SP_source_3d &source)
{
  Require(source);

  // Build the fission density for keff = 1.1
  source->build_fission_density(1.1);

  // Get the density.
  const FissionSource<DD_3D_Equations,_3D>::FissionDensityField &f =
      source->get_fission_density();
  FissionSource<DD_3D_Equations,_3D>::SP_mesh mesh = source->mesh();
  // Verify length of density vector
  TEST(f.size() == mesh->number_cells());

  // Check fission density.
  const Mesh<_3D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(f[i], f_ref[matid[i]]));

  cout << " ... passed 3-d FissionSource build_fission_density test. " << endl;
  return 0;
}

int test3d_build_fission_source(SP_source_3d &source)
{
  Require(source);

  // Build the fission density for group 1 (indexed as 0!!)
  source->build_fission_source(0);

  // Get the source.
  const FissionSource<DD_3D_Equations,_3D>::FissionSourceField &f =
      source->get_fission_source();
  FissionSource<DD_3D_Equations,_3D>::SP_mesh mesh = source->mesh();
  // Verify length of density vector
  TEST(f.size() == mesh->number_cells());

  // Set reference factor.  Chi is 1, so factor is just 1/(4pi*keff) for 2d.
  double factor = inv_four_pi / 1.1;

  // Check factor.
  TEST(source->factor() == factor);

  // Check fission source vector directly.
  const Mesh<_3D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(f[i], factor * f_ref[matid[i]]));

  cout << " ... passed 3-d FissionSource build_fission_source test. " << endl;
  return 0;
}

int test3d_verify_source(SP_source_3d &source)
{
  Require(source);

  // Build the fission density for group 1 (indexed as 0!!)
  source->build_fission_source(0);

  // Get the mesh.
  FissionSource<DD_3D_Equations,_3D>::SP_mesh mesh = source->mesh();

  // Set reference factor.  Chi is 1, so factor is just 1/(4pi*keff) for 3d.
  double factor = inv_four_pi / 1.1;

  // Check factor.
  TEST(source->factor() == factor);

  // Check fission source using q_f interface.
  const Mesh<_3D>::Vec_Id matid = mesh->get_cell_id(MATERIAL_ID);
  for (int i = 0; i < mesh->number_cells(); i++)
    TEST(soft_equiv(source->q_f(i, 0, 0), factor * f_ref[matid[i]]));

  cout << " ... passed 3-d FissionSource verify_source test. " << endl;
  return 0;
}

//---------------------------------------------------------------------------//
//              end of testFissionSource.cc
//---------------------------------------------------------------------------//
